function XglomRunPipeline(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%
%   - Michael Eager,   (michael.eager@monash.edu) 
%   - (c) 2012, Monash Biomedical Imaging, Monash University, Australia

%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

display 'Running "Xglom Pipeline version 1.6" ... '

debug=1;
splitkidney=0;

% MR pipeline variables
Xglom_MR_pipeline.version=handles.Xglom_version;
Xglom_MR_pipeline.ProcessingDate=datestr(now);
Xglom_MR_pipeline.CorrMask_threshold=1.75;
Xglom_MR_pipeline.CorrMask_sigma=3.0/60.0;
Xglom_MR_pipeline.RegularizeMRI_factor=1;  %default
Xglom_MR_pipeline.RegularizeMRI_filter_sd=2.0; %default




tstart=tic;
handles.Img = GetDICOMImage(handles.DatasetID);
addpath(fullfile(pwd,'plugins','progressbar'));
gui_active(1);      % will add an abort button
hDummy=figure;hProgress=figure;
hProgress_0 = hProgress;hProgress  = progressbar([],0,' Phase 1 of 2: Image Cleaning Procedure (8 steps)' );
ShowStatus('Xglom Pipeline : Phase 1 of 2: Image Cleaning Procedure (8 steps)',handles);incr_progress = double(1)/double(60); %estimate 240 secs for imfill, 50 reps for bwconncomp
midaxis = floor(size(handles.Img,1)/2);



%disable original hProgress figure
close(hProgress_0);
close(hDummy)


%% Phase 1a: Grab grayscale mask
ShowStatus('Phase 1a: Grab bias mask',handles);
% Use default args for RegularizeMRI.  To set slice range  499 to 504, with
% factor = 1 use mask=RegularizeMRI(img,499,504,1) 
mask = RegularizeMRI(handles.Img);
if debug==1
    %fig=figure;imshow(mask);pause; close(fig);
    imagesc(mask,'Parent',handles.pipelineaxes)
end
%% Phase 1b: Apply grayscale mask to volume
hProgress = progressbar( hProgress, incr_progress, 'Phase 1: Step 2 of 8 - Grayscale mask');
img_clean = MaskVolume(img,mask);
if debug==1
    fig=figure;imagesc(squeeze(img_clean(256,:,:)));colormap(flipud(gray));pause; close(fig);
end

%% Phase 1c: Apply Laplacian of Gaussian filter to cleaned image
hProgress = progressbar( hProgress,incr_progress, 'Phase 1: Step 3 of 8 - L-o-G filter');
h = fspecial3('log',5);
Ilog = imfilter(img_clean,h,'replicate');
%imshow(squeeze(Ilog(256,:,:)))
if debug==1
    fig=figure;imagesc(squeeze(Ilog(256,:,:)));pause; close(fig);
end


%% Phase 1d:  Smooth out Ilog greater than zero with inhomogenityCorrection
hProgress = progressbar( hProgress,incr_progress, 'Phase 1:  Step 3 of 8 - Smooth out L-o-G');
tic;corrMask=inhomogenityCorrection((Ilog.^2).*(Ilog>0),3./60);toc
%imshow(squeeze(corrMask(256,:,:)))
if debug ==1
fig=figure;imagesc(squeeze(corrMask(256,:,:)));pause; close(fig);
end

%% Phase 1e: Set a mask in smoothed-out Ilog and find connected components
hProgress = progressbar( hProgress,incr_progress, 'Phase 1:  Step 4 of 8 - Find glomerular regions');
CC = bwconncomp(corrMask>1.75);  % threshold is selected from histogram
numPixels = cellfun(@numel,CC.PixelIdxList);
[x idx]=sort(numPixels,'descend');
idx_lim=1;

if  sum(x(1:idx_lim)) < 15000000   
% 15000000 = 0.111*512^3
     disp('Binary labels not big enough to encompass two kidney glomerular regions, increasing regions.')
    idx_lim=idx_lim+1;
end

%% Phase 1f:  Create binary mask
hProgress = progressbar( hProgress, incr_progress, 'Phase 1:  Step 5 of 8 - Create binary mask');
BW = ismember(labelmatrix(CC), idx(1:idx_lim));
%  9.016970 seconds.
%imshow(squeeze(BW(256,:,:)))
if debug==1
fig=figure;imagesc(squeeze(BW(256,:,:)));pause; close(fig);
end




%% Phase 1g: Image filling
hProgress = progressbar( hProgress, incr_progress, 'Phase 1:  Step 6 of 8 - Flood-filling');
telapsed = toc(tstart)
imf=single(imfill(img_clean.*BW,'holes')); %default:conn=6
depth=imf - (img_clean.*BW);

%imshow(squeeze(depth(256,:,:)))
if debug==1
fig=figure;imagesc(squeeze(depth(256,:,:)));pause; close(fig);
end




%% Phase 1h: Normalise the depth image and save processed data
hProgress = progressbar( hProgress, incr_progress, 'Phase 1:  Step 7 of 8 - Normalising depth image');
depth = NormaliseImage(depth);

system(['[ -f ./Cache/' handles.DatasetID '/depth_mask_01.mat ] && chmod ug+rw ./Cache/' handles.DatasetID '/depth_mask_01.mat']); 
save(['./Cache/' handles.DatasetID '/depth_mask_01.mat'],'depth','mask');
system(['chmod ug+rw ./Cache/' handles.DatasetID '/depth_mask_01.mat']); 


%% Phase 2: Thresholding
hProgress = progressbar( hProgress, incr_progress, 'Phase 2: Performing Threshold Statistics (50 steps)');

telapsed = toc(tstart)
display 'Running Statistical tests ...'

ThresholdRange= 0.01:0.01:0.5;
stored_hist=zeros(50,1000);
stored_diams = zeros(50,141);
counter=1;
for counter=1:length(ThresholdRange) % Loop over different Thresholds
  %% Thresholding
  Holes=depth>ThresholdRange(counter);

  %% Discard single voxels
  % ConArray=conndef(3,'minimal');
  % ConArray(2,2,2)=0;
  % Holes=imdilate(Holes,ConArray).*Holes;
  % Holes=single(Holes);

  %% Statistics, Label Glomerular Regions
  CC = bwconncomp(Holes,6);
  if counter == 20 
    labeledgloms =  labelmatrix(CC);
  end
  
  % calculate area histogram
  areas = cellfun(@numel,CC.PixelIdxList);
  area_max = 1000; %max(areas);
  area_range=(1:area_max);
  [area_hist, area_axis]=hist(areas,area_range);
  stored_hist(counter,:) = area_hist;


  % Total number of detected objects in area_range
  ['sum(area_hist(1:100)):  ', num2str(sum(area_hist(1:100)))]
  
  % Diameter histogram from area data assuming spherical glomeruli
  [diams x]=hist((6/pi*areas(areas>2)).^(1/3),1:0.1:15); diams(end)=0;
  stored_diams(counter,:) = diams;
  

  if counter == 50
    hProgress = progressbar( hProgress, incr_progress, ['Phase 2: Completed Xglom pipeline. Saving Data.']);  
  else
    hProgress = progressbar( hProgress,incr_progress , ['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)']);  
  end
end
telapsed = toc(tstart)

%% Write output data to file - ensure writability by users and group
system(['[ -f ./Cache/' handles.DatasetID '/output_01.mat ] && chmod ug+rw ./Cache/' handles.DatasetID '/output_01.mat']); 
save(['./Cache/' handles.DatasetID '/output_01.mat'],'stored_hist', 'stored_diams', 'labeledgloms')
system(['chmod ug+rw ./Cache/' handles.DatasetID '/output_01.mat']); 


hProgress = progressbar( hProgress, -1 , ['Image Processing Complete. Processing output plots.']); 

PlotXGlomOutput(stored_hist,stored_diams,handles.DatasetID);
handles.stored_hist = stored_hist;
handles.stored_diams = stored_diams;

